#ifndef UTILS_H
#define UTILS_H
#include <complex>
#include <x86intrin.h>
#include <chrono>
// #include <mkl.h>
// #include <omp.h>
typedef std::chrono::high_resolution_clock Clock;

// #define THE_THREAD_NUMS 8
using namespace std;
// #define NO_MPI_IO
extern uint64_t count_mpi_time01;
extern uint64_t count_mpi_time02;
extern uint64_t count_mpi_time03;
extern int onlyonce_01;
extern int onlyonce_02;
extern int onlyonce_03;
extern int global_rank;

inline void cal_x_f(complex<double> *srcO, complex<double> *U, complex<double> *destE, double flag,
                    complex<double> I)
{
    std::complex<double> tmp;
    const double half = 0.5;
    for (int c1 = 0; c1 < 3; c1++) {
        for (int c2 = 0; c2 < 3; c2++) {
            {
                tmp = -(srcO[0 * 3 + c2] - flag * I * srcO[3 * 3 + c2]) * half * U[c1 * 3 + c2];
                destE[0 * 3 + c1] += tmp;
                destE[3 * 3 + c1] += flag * (I * tmp);
                tmp = -(srcO[1 * 3 + c2] - flag * I * srcO[2 * 3 + c2]) * half * U[c1 * 3 + c2];
                destE[1 * 3 + c1] += tmp;
                destE[2 * 3 + c1] += flag * (I * tmp);
            }
        }
    }
}

inline void cal_x_b(complex<double> *srcO, complex<double> *U, complex<double> *destE, double flag,
                    complex<double> I)
{
    std::complex<double> tmp;
    const double half = 0.5;
    for (int c1 = 0; c1 < 3; c1++) {
        for (int c2 = 0; c2 < 3; c2++) {
            tmp = -(srcO[0 * 3 + c2] + flag * I * srcO[3 * 3 + c2]) * half * conj(U[c2 * 3 + c1]);
            destE[0 * 3 + c1] += tmp;
            destE[3 * 3 + c1] += flag * (-I * tmp);

            tmp = -(srcO[1 * 3 + c2] + flag * I * srcO[2 * 3 + c2]) * half * conj(U[c2 * 3 + c1]);
            destE[1 * 3 + c1] += tmp;
            destE[2 * 3 + c1] += flag * (-I * tmp);
        }
    }
}

// inline void cal_y_f(complex<double> *srcO, complex<double> *U, complex<double> *destE, double flag,
//                     complex<double> I)
// {
//     std::complex<double> tmp;
//     const double half = 0.5;
//     for (int c1 = 0; c1 < 3; c1++) {
//         for (int c2 = 0; c2 < 3; c2++) {
//             tmp = -(srcO[0 * 3 + c2] + flag * srcO[3 * 3 + c2]) * half * U[c1 * 3 + c2];
//             destE[0 * 3 + c1] += tmp;
//             destE[3 * 3 + c1] += flag * (tmp);
//             tmp = -(srcO[1 * 3 + c2] - flag * srcO[2 * 3 + c2]) * half * U[c1 * 3 + c2];
//             destE[1 * 3 + c1] += tmp;
//             destE[2 * 3 + c1] -= flag * (tmp);
//         }
//     }
// }

inline void cal_y_f(complex<double> *srcO, complex<double> *U, complex<double> *destE, double flag,
                    complex<double> I)
{
    std::complex<double> tmp;
    std::complex<double> ss_0[3];
    std::complex<double> ss_1[3];
    for (int c2 = 0; c2 < 3; c2++) {
        ss_0[c2] = srcO[0 * 3 + c2] + flag * srcO[3 * 3 + c2];
        ss_1[c2] = srcO[1 * 3 + c2] - flag * srcO[2 * 3 + c2];
    }
    const double half = 0.5;
    for (int c1 = 0; c1 < 3; c1++) {
        for (int c2 = 0; c2 < 3; c2++) {
            tmp = -(ss_0[c2]) * half * U[c1 * 3 + c2];
            destE[0 * 3 + c1] += tmp;
            destE[3 * 3 + c1] += flag * (tmp);
            tmp = -(ss_1[c2]) * half * U[c1 * 3 + c2];
            destE[1 * 3 + c1] += tmp;
            destE[2 * 3 + c1] -= flag * (tmp);
        }
    }
}

inline void cal_y_b(complex<double> *srcO, complex<double> *U, complex<double> *destE, double flag,
                    complex<double> I)
{
    std::complex<double> tmp;
    const double half = 0.5;
    for (int c1 = 0; c1 < 3; c1++) {
        for (int c2 = 0; c2 < 3; c2++) {
            tmp = -(srcO[0 * 3 + c2] - flag * srcO[3 * 3 + c2]) * half * conj(U[c2 * 3 + c1]);
            destE[0 * 3 + c1] += tmp;
            destE[3 * 3 + c1] -= flag * (tmp);
            tmp = -(srcO[1 * 3 + c2] + flag * srcO[2 * 3 + c2]) * half * conj(U[c2 * 3 + c1]);
            destE[1 * 3 + c1] += tmp;
            destE[2 * 3 + c1] += flag * (tmp);
        }
    }
}

inline void cal_z_f(complex<double> *srcO, complex<double> *AE, complex<double> *destE, double flag,
                    complex<double> I)
{
    std::complex<double> tmp;
    const double half = 0.5;
    for (int c1 = 0; c1 < 3; c1++) {
        for (int c2 = 0; c2 < 3; c2++) {
            tmp = -(srcO[0 * 3 + c2] - flag * I * srcO[2 * 3 + c2]) * half * AE[c1 * 3 + c2];
            destE[0 * 3 + c1] += tmp;
            destE[2 * 3 + c1] += flag * (I * tmp);
            tmp = -(srcO[1 * 3 + c2] + flag * I * srcO[3 * 3 + c2]) * half * AE[c1 * 3 + c2];
            destE[1 * 3 + c1] += tmp;
            destE[3 * 3 + c1] += flag * (-I * tmp);
        }
    }
}

inline void cal_z_b(complex<double> *srcO, complex<double> *AO, complex<double> *destE, double flag,
                    complex<double> I)
{
    std::complex<double> tmp;
    const double half = 0.5;
    for (int c1 = 0; c1 < 3; c1++) {
        for (int c2 = 0; c2 < 3; c2++) {
            tmp = -(srcO[0 * 3 + c2] + flag * I * srcO[2 * 3 + c2]) * half * conj(AO[c2 * 3 + c1]);
            destE[0 * 3 + c1] += tmp;
            destE[2 * 3 + c1] += flag * (-I * tmp);
            tmp = -(srcO[1 * 3 + c2] - flag * I * srcO[3 * 3 + c2]) * half * conj(AO[c2 * 3 + c1]);
            destE[1 * 3 + c1] += tmp;
            destE[3 * 3 + c1] += flag * (I * tmp);
        }
    }
}
inline void cal_t_f(complex<double> *srcO, complex<double> *AE, complex<double> *destE, double flag,
                    complex<double> I)
{
    std::complex<double> tmp;
    const double half = 0.5;
    for (int c1 = 0; c1 < 3; c1++) {
        for (int c2 = 0; c2 < 3; c2++) {
            tmp = -(srcO[0 * 3 + c2] - flag * srcO[2 * 3 + c2]) * half * AE[c1 * 3 + c2];
            destE[0 * 3 + c1] += tmp;
            destE[2 * 3 + c1] -= flag * (tmp);
            tmp = -(srcO[1 * 3 + c2] - flag * srcO[3 * 3 + c2]) * half * AE[c1 * 3 + c2];
            destE[1 * 3 + c1] += tmp;
            destE[3 * 3 + c1] -= flag * (tmp);
        }
    }
}
inline void cal_t_b(complex<double> *srcO, complex<double> *AO, complex<double> *destE, double flag,
                    complex<double> I)
{
    std::complex<double> tmp;
    const double half = 0.5;
    for (int c1 = 0; c1 < 3; c1++) {
        for (int c2 = 0; c2 < 3; c2++) {
            tmp = -(srcO[0 * 3 + c2] + flag * srcO[2 * 3 + c2]) * half * conj(AO[c2 * 3 + c1]);
            destE[0 * 3 + c1] += tmp;
            destE[2 * 3 + c1] += flag * (tmp);
            tmp = -(srcO[1 * 3 + c2] + flag * srcO[3 * 3 + c2]) * half * conj(AO[c2 * 3 + c1]);
            destE[1 * 3 + c1] += tmp;
            destE[3 * 3 + c1] += flag * (tmp);
        }
    }
}

#endif
